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ABSTRACT 

Stars near the Galactic center reach a few percent of light speed during peri- 
center passage, which makes post-Newtonian effects potentially detectable. We 
formulate the orbit equations in Hamiltonian form such that the 0(f^/c^) and 
0{v^/c^) post-Newtonian effects of the Kerr metric appear as a simple general- 
ization of the Kepler problem. A related perturbative Hamiltonian applies to 
photon paths. We then derive a symplectic integrator with adaptive time-steps, 
for fast and accurate numerical calculation of post-Newtonian effects. Using this 
integrator, we explore relativistic effects. Taking the star S2 as an example, we 
find that general relativity would contribute tenths of mas in astrometry and 
tens of km s~^ in kinematics. (For eventual comparison with observations, red- 
shift and time-delay contributions from the gravitational field on light paths will 
need to be calculated, but we do attempt these in the present paper.) The con- 
tribution from stars, gas, and dark matter in the Galactic center region is still 
poorly constrained observationally, but current models suggest that the resulting 
Newtonian perturbation on the orbits could plausibly be of the same order as 
the relativistic effects for stars with semi-major axes > 0.01 pc (or 250 mas). 
Nevertheless, the known and distinctive time dependence of the relativistic per- 
turbations may make it possible to disentangle and extract both effects from 
observations. 
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Introduction 



The Galactic-center stars are a population of fast-moving stars in highly eccentric nearly 
Keplerian orbits around a compa ct mass — presumably a massive black hole (hencefort h 
MBH) — of mass ~ 4 x IO^Mq Jchez et al.lboosl : Ichez et alj[2008l : Icillessen et aDboOoh . 



Pericenter distance for some of these stars are inferred to be < 10^ of the gravitational radius 
(~ 3 X 10^ in the case of S2) implying pericenter velocities of a few percent of light. 

The high pe ricenter velo c ities i nspir e a search for general rela tivistic perturbations to 
Keplerian orbits. IJaroszynskil (119981 ) and [Fragile fc Mathewd (120001 ) suggested that the pre- 
cession of the pericenter, which is an effect of 0{v'^) where v is the pericenter velocity in light 
units, would be observable. Interfe rometric instruments currently under development (see, 
for example, lEisenhauer et al.ll2009l ) make this possibility more likely. If stars even further in 
are discovered, then 0(v^) could become detectable from astrometry, leading to tests of no- 
hair theorems ( Will 20081 ). O ther consequences o f general relativity may become accessible 
to spectroscopic observations. IZucker et al.l (120061 ) note that the O(v^) ef f ect of gravitational 
redshift may be within the reach of current instruments. iKannan fc Sahal (120091 ) suggest that 
the next generation of spectrographs may pick up even the O(w^) signature of frame dragging 
for the known Galactic-center stars. 

To study the relativistic perturbations in detail, good methods for computing them 
numerically are highly desirable. In this paper we provide such a method. 



Since geodesies in a Kerr metric are integrable (see, for example, chapter 7 of I Chandr asekhar 



19831 ). our problem may at first seem a trivial one. But the known integrability only puts 
the solution in terms of quadratures, it does not p rovid e solutions in terms of elementary 
functions. The explici t forms given by iKraniotid (120071 ) for polar orbits, and the numer- 
ical methods given by iDexter fc Agoll (120091 ) for computing null geodesies give an idea of 
the complexity of the formally exact solutions. Moreover, even if one had the exact Kerr 
geodesies, one would need to perturb them for the Galactic-center stars, because of masses 
other than the black hole. 

In this paper we take a different route. We go back to the Kerr metric and derive a 
post-Newtonian Hamiltonian, where relativistic effects appear as perturbations to Keplerian 
orbits and to which Newtonian perturbations due to other masses can be trivially added. 
Then we consider numerical algorithms designed for comets and other highly eccentric orbits 
in the solar system, and generalize them to work for post-Newtonian perturbations as well. 
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2. A post-Newtonian Hamiltonian 



Geodesic equations in a given metric can be described in terms of a Lagrangian 



2 dtif 



dr dr ' 



(i; 



where r is the afhne parameter which, conveniently for us, can be identified with the star's 
proper time. Less common is the equivalent Hamiltonian form 



(2) 



In fact H and L are equal, and conserved along geodesies. The above duality is simply a 
consequence of being quadratic. 

We want to derive a post-Newtonian Hamiltonian for a test particle in a Kerr metric. 
To do this, we write the metric in Boyer-Lindquist coordinates, adopting (— , +, +, +) for 
the signature of the metric and choosing units in which G = c = 1, as follows: 
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where the dimensionless spin parameter s G [0, 1]. The contravariant components are, as is 
not too difficult to verify by inverting the matrix, as follows: 
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Expanding to second order in /x, a, and n we have 



\ 



-(1 + 2/X + V 



2/i(j 



1 - 2/i + CT^ 
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sin^ ^ ^ 



Now we proceed to consider the Hamiltonian. For simphcity, we put M = 1, which just 
means we are measuring r in units of the gravitational radius, GM/c^. To consider dynamics 
at large-r we replace 

r 



The Hamiltonian will now have two regimes. 

In the low- velocity regime we require the velocity terms to be 0(e). This is achieved 
with the substitutions 



Pr 



epr 



which gives 

2 



2 V 2 2r2 + 2r2 sin^ 6 
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(9) 



-e^ + 0(e^). (10) 



At zeroth order, a test particle just stays still. At ©(e^) it follows Newtonian dynamics. 
Relativistic effects appear at 0{e'^), while frame-dragging appears at O(e^). But note that 
the kinematic effects themselves are at one order of e lower: thus Newtonian velocities 
are 0(e), Schwarzschild perturbations to the velocities are O(e^), while the frame-dragging 
perturbation to velocity is O(e^). Gravitational radiation is a higher-order effect which we 
disregard here, and in any case the timescale for orbital decay of the S2 star due to radiation 
reaction is longer than a Hubble time. 

In the light- velocity regime we require the velocity terms to be 0(1). Thus we replace 



which gives 



H = Pr^ Pe'^ 

2 2 2r2 2r2sin2e 



2pt^ _ sin^ e 



2 S^COS^6' 2 

Pr + + 



^1 e^ + 0(e^) 



(12) 



^Here e is just a label for keeping track of orders. Numerically e = 1. 
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At zeroth order, null geodesies just move in straight lines. The leading perturbation is at 

The Hamiltonian flTUl) is the approximate Hamiltonian we will use, but we can simplify 
its form with some variable changes. First, we set = —1, which we are free to do since 
the Hamiltonian is autonomous. This merely sets units for the affine parameter such that 
dt/dr = —1 in the large-r limit, and has no physical significance. Then, we change from 
r, 6, (f) to X, y, z. Completing the canonical transformation, we have 

Pr = , P0 = (X X p)^, (13) 

r 

and hence 

+ ^ + ^ = p'- (") 

sm y 

The post-Newtonian Hamiltonian ffTOl) then becomes 

H = Hkcp + Hs + Hut, (15) 

where 

n 2 (x-p)^ 

= --^ (16) 



There is a separate equation for t 



S ■ X X p 



t=i + H + l_H£|*. (17) 



3. An adaptive-timestep symplectic integrator 

The post-Newtonian orbit equations can be integrated numerically by any general- 
purpose method for ordinary differential equations. Another option is to use an iV-body sim- 



ulatio n code for dense stellar systems, with post-Newtonian terms added (IMikkola fc Merritt 



20081 ). A computationally more efficient strategy, however, would be an integration algorithm 
that takes advantage of the above formulation of relativistic effects as small perturbation 
to a Kepler Hamiltonian. We now design such an integration algorithm, based on recent 
work on cometary orbits, which are also highly eccentric orbits that experience interesting 
perturbations around pericenter passage. 
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When integrating Hamiltonian systems numerically, it is a common practice to impose 
the condition that the numerical solution has (to machine precision) the symmetry properties 
of Hamiltonian flow. Integration algorithms with this property are known as symplectic 
integrators. A simple but important example is generalized leapfrog. Suppose we have a 
Hamiltonian that is the sum of two parts H = Ha + Hb where Ha and Hb are individually 
easy or trivial to integrate. Generalized leapfrog evolves under H for a time step At as 
follows. 

1. Evolve under Ha for time |Ar. 

2. Evolve under Hb for time At. 

3. Reiterate step 1. 



If if = + ^(x) the above bec omes the classical leapfrog integrator. It turns out (see, 



for example, iForest fc Ruth! 119901 ) that generalized leapfrog amounts to evolving under a 
"surrogate" Hamiltonian 

Ha + Hb+ [U^a, Hb}, Hb} + 1,{Ha, Hb}, Ha})At^ + O(Ar^). 

[ Poisson brackets amount to a Hamiltonian expression for the erro 



The nested Poisson brackets amount to a Hamiltonian expr ession for the error, which is 
mani f estly second order. High er-order extensions are possible (IForest fc Ruthlll990l : lYoshida 
1990l : iLaskar &: Robutell l200ll ) but in practice second-order is the most used. If one of Ha 
or Hb is much smal l er tha n the other, the erro r Hamiltonian w i ll be correspondingly small. 
Wisdom fc HolmanI (119911 ) and independently iKinoshita et al.l (jl99ll ) proposed integrators 
for planetary orbits where Ha is the integrable Kepler Hamiltonian and Hb encapsulates 
the perturbations. For planetary orbits, the low order of generalized leapfrog becomes an 
advantage in that the steps can be made comparatively large (~ 10 steps per orbit) and 
still provide high accuracy. Furt her refinements are pos sible, such as perturbative pre- 
processing of t he initial conditions (ISaha &: Tremaindll992l ) or perturbative post-processing 
of the results (IWisdom et al.lll996l ). but the original Wisdom-Holman scheme is the most 
common choice for long-term solar-system orbit integrations. 

Another ki nd of symp lectic integrator, also second order, is the implicit midpoint 
method, due to iFengI (119861 ). This can be written as a simple discretization of Hamilton's 
equations 

'dH\ ^ ^ fdH\ 



Ax 



with the derivatives are evaluated at the midpoint 

(p + iAp,x+iAx) 



(19) 



(20) 
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It is essential for Ax, Ap to be consistent between f|T9|) and fl20l) to high accuracy (preferably 
machine precision), otherwise the symplectic property is lost. Hence, the implicit midpoint 
method requires iteration. But it requires no splitting of the Hamiltonian, and hence is very 
useful when generalized leapfrog is inapplicable. A surrogate Hamiltonian for the implicit 
midpoint integrator is derived in ISaha et al.l (119971 ). 



Combining the two above ingredients, a possible integration method for the Galactic- 
center stars would be a generalized leapfrog with Ha = -ffKep and Hb = + i^LT + Vcai, 
with an exact Kepler solution used for the former, and implicit midpoint used for the latter. 
In fact, neither of the separate integrations under Ha and Hb needs to be exact. As long as 
they are symplectic and second-order, the surrogate Hamiltonian f|T8l) will apply. 

For low eccentricites, an integrator as above would be very efficient. When we consider 
Galactic-center stars (or comets) however, we run into the major limitation of generalized 
leapfrog: the stepsize At must remain fixed, otherwise the integrator is no longer symplectic. 
Yet for e ~ 0.9, the very small At needed at pericenter becomes hopelessly expensive if used 
throughout an orbit. Adaptive time-stepping is needed. 

The key to adaptive time-stepping is to transform from r to a new independent vari- 
able (say s) which somehow implements the des ired stret c hing and shrinking of the step- 



size without breaking the symplectic property. iMikkolal (119971 ) provided the first exam- 
ple, showing how a simple modification of the Wisdom-Holna a n alg orithm effective cre- 
ates the variable s with dT = rds. Then iPreto &: Tremaind (119991 ) and independently 



Mikkola fc Tanikawal (119991 ) formulated, for a Hamiltonian ^p^ + V(x), a time transforma- 
tion with dT = —ds/V{x^). As a by-product, this work produced a leapfrog for H^ep that 
is symplectic and in fact recovers the exact answer except for a very small pha s e erro r, but 
is computationally much simpler than the exact solution. iMikkola fc AarsethI (120021 ) then 
generalized these ideas to make the time-transformation completely a daptive, though not in 
a Hamiltonian formulation, so it was not manifestly symplectic. Later, lEmeFyanenkd (120071 ) 
supplied an elegant Hamiltonian derivation of adaptive stepsize. 

Based on al l the above, w e now develop our integrations algorithm. The derivation 
basically follows lEmePyanenkd (120071 ) but is written with a view to application to post- 
Newtonian orbits. 

Let us enhance the phase space (x, p). We now treat r as an additional coordinate, 
with a new variable $ being its conjugate momentum, and the new independent variable s. 
In this enhanced phase space, consider the Hamiltonian 



F(x,r,p,$) 



H -k 



+ In 



¥?(r) 



(21) 
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where (p{t) is a known function and is a constant. Hamilton's equations for F are 

dx _ IdH dp _ IdH 
ds ^ dp ds $ (9x ' 

together with 



(22) 



:r = :r^^^ (23) 

ds dr 

and 

dr H-k 1 

Moreover, F will be conserved. 

Suppose at some s we have H = k. Combining the two previous equations, we infer 
that at this point 

In $ = ^ In (25) 
dr dr 

Thus, in the neighborhood of this point, will be constant, and hence H — k will remain 
zero. In other words, if 

holds initially, it will continue to hold. The interpretation is the original Hamiltonian equa- 
tions with a rescaled variables ds = ^ dr. 

Now the cunning part: we consider ^plr) as fl'x, p) where x, p are functions of r, and 
replace dip/dr by the convective derivative. We have: 

d^ f dy. d dp d \ , , ^ , , 

With the above ingredients in hand, we proceed to write a leapfrog for F in the variable 

s: 



1. Holding X, p constant, advance $ using (!27|) over As = |. This amounts to evolution 
under a Hamiltonian — Imp. 

2. Holding $ constant, evolve x, p under H/^ for As = 1, and advance r by 1. This 
amounts to evolution under a Hamiltonian {H — A;)/<l> + ln$. 

3. Reiterate step 1. 



- 9 - 



In the present work, step 2 above consists of a generalized leapfrog with the main Hamiltonian 
split into Keplerian, post-Newtonian, and external potential parts. We also have two minor 
simplifications: (i) we assume ip depends only on x, and (ii) we disregard the evolution of r. 
The simplified treatment of time is fine for our test particle integrations, but for an A^-body 
formulation, the time equation needs to be treated with special care, and the form dr = ds/^, 
associated with the new pair of conjugate variables (r, $) presents several advantages. 

The algorithm is as follows. 



1. Advance $ by ^{dH/dp) ■ V(ln(/?). 

2. E volve X, p under Hxep for At = 1/(2$), using the algorithmic regularization scheme 
of IPreto fc Tremainel fll999l ). 

3. Evolve X, p under + i^LT + Vcai for = V*^; using the implicit midpoint method. 
For the small post-Newtonian contributions of interest for this paper, this implicit 
scheme converges to machine precision in two or there iterations. 

4. Reiterate step 2. 

5. Reiterate step 1. 



The following comments about the integrator are worth making. First, with the simplest 
choice (/? oc 1/r, we obtain an adaptive, symplectic integrator to integrate eccentric weakly- 
perturbed Keplerian orbits — including non-separable post-Newtonian p erturbations. This 
i ntegr ator is easy to implement, and it is free from the instability found iRauch fc Holman 
( I1999I ) in fixed-stepsize integration of highly eccentric orbits. Second, with the apparently 
innocuous modification of the time transformation to dr = ds/^, the Keplerian part be- 
comes trivial to integrate, since $ is kept constant while the (cartesian) coordinates are 
advanced. This is very advantageous for a A^-body implementation with individual time 
steps, where particle synchronization requires numerous Kepler drifts. Third, being a sym- 
plectic integrator, errors in the longitudes grow only linearly with time and, for spherical 
perturbations, angular momentum is conserved to machine precision. This is clearly a most 
desirable property — in particular, in dynamical problems for which long term integrations 
with accurat e tracking of all phase angles are required, e. g., resonant relaxation and re- 



lated effects (iGiirkan fc Hopmanl [2OO7I : Madigan et allbood : IPerets et al.lbood l Fourth 



as 



pointed out already by Emel'yanenkol ( 2007 ). the freedom to choose the functional form of 
if provides an additional degree of freedom that can be explored in order to resolve close 
encounters without breaking symplecticity. Fifth, this integrator and the post-Newtonian 
approximation presented in this work are very well suited to be included in orbital fitting 
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routines of the Galactic-center stars (jWeinberg et al.ll2005l : iGillessen et al.ll2009l ). Sixth, this 
symplectic integrator is also ideally suited for implementation in view of A^-body modelling 
of gravitational wave sources, e.g. extreme mass ratio inspirals, which are of great interest 
for LISA (Laser Interferometer Space Antenna). 



Numerical results 



As an illustration, we consider perturbations of the orbit with Keplerian elements 
a = 2.4 X 10^ e = 0.88 / = 135.25° VL = -134.71° uj = 63.56° (28) 



which are approximately the measured values for S2 (IGillessen et al.l 120091 ). Note that a is 
in units of the gravitational radius GM/c^ ~ 5 x 10^ km or ~ 4/ias on the sky. We assume 
the spin is unity (maximal) and directed along +z. 

We begin by verifying that the integrator is indeed second order and that the phase 
error grows linearly with time. From Figure [1] we see that the maximum error in H is 
proportional to $"^(0), where 1/<I'(0) is the initial stepsize, and that the error in pericenter 
angle is proportional to the number of orbits [i.e. it is linear rather quadratic). In order to 
evaluate this error we ado pt the expressions for the shift of pe ricenter due to a Schwarzschild 



black hole up t o 2 order (iWeinberg 



19721 : iNucita et al.ll2007l ) and the l"** order contribution 



from the spin (IKannan fc Sahall2009l ): 



67r 



a 1 



+ 



37r(18 + e2l 



2a2(i-e2)2 



(29) 



Aa; 



127rs cos/ 



^3/2(1 



^2^3/2 ■ 



The observable orbit, that is to say, the sky position and redshift as a function of 
observer time, depends also on the light path from the star to the observer. There are two 
types of effects. 



1. The sky position of the star will be slightly shifted by gravitational lensing. The 
maximum lensing displacement is the Einstein radius Re. Since for Galactic- Center 
stars, the lens-source distance Dls is much smaller than the observer-source and lens- 
source distances, the usual expression for Re simplifies to 

Re = (30) 

in units of the gravitational radius, and D^s will be of order a. 
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2. Then there is the redshift, and the R0mer effect, which can be considered as the redshift 
integrated along the orbit. The 0(e) contribution to the redshift is classical. Special 
relativity and the time part of the metric both contribute to the redshift at O(e^), 
while space curvature contributes at O(e^). 



We do not include all these effects in this paper, leaving it for future work, because at this 
stage our aim is to gain some insight into the size of relativist ic effects rather than calculate 



observables precisely. However, following IZucker et al.l (120061 ) . we include the O(f^) effects 



due to gravitational time dilation when we estimate the perturbations on the radial velocities. 

Having chosen the orbital elements, we start the star at apocenter and integrate the 
orbit under the Hamiltonian f[T^ to the next apocenter passage. Figure [2] shows a, u;, and Q 
along the orbit, as a function of the mean anomaly, which is a surrogate for time. The range 
— 180 to 180 in the mean anomaly corresponds to a Keplerian orbit going from apocenter 
to apocenter. In this paper we always compare perturbed and unperturbed orbits for the 
same value of mean anomaly, not necessarily the same value of time. If time is taken as the 
independent variable, an artificial secular drift would appear, because the perturbed orbit 
has a slightly different orbital period. The other two Keplerian elements are not shown here, 
because / and a(l — e^) are constant as a consequence of the conservation of total angular 
momentum, which the symplectic integrator reproduces to roundoff error. 



From Figure [21 one easily verifies the well-known leading-order expressions 

ail — e 

for the pericenter precession and 



= 77^— IT (31) 



. ^ Stt cos / , , 

An = ^ 32 

for Lense-Thirring node precession. The maximum change in a can be estimated as follows. 
We start by noting from the inspection of f|T6l) that the leading post-Newtonian perturbation 

2 

^H=^ (33) 
at both pericenter and apocenter. Since at these points r = a(l ± e) we have 



e 

a'Hl -e\ 

2\ 



Ai/pe. - = -^7—-^ (34) 



For the unperturbed Hamiltonian H = — l/(2a) hence AH ~ Aa/(2a ) from which it follows 
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Using similar arguments, or dimensional analysis, one easily derives that the effects of 
a perturbation 

~ r"" (36) 

scale as follows. 

astrometric A/, A6 ~ a^~" 

kinematic Av ~ a^""' (37) 

precession Acj, etc ~ a^^"' 

Post-Newtonian effects have n = 2 (Schwarzschild) or = 3 (Lense-Thirring). Thus, the 
astrometric effect of Schwarzschild perturbations is independent of a to leading order, while 
other post-Newtonian effects get stronger as orbits get smaller. Relativistic prograde pre- 
cession would seem easier to measure at larger distances, but that is not the case since the 
orbital period increases as a?/"^ . On the other hand. Galactic perturbations from other masses 
than the black hole would have n < 1, and hence get weaker as orbits get smaller. 

The form and strength of Newtonian Galactic perturbations, due to other stars, gas, 
and dark matter, are still poorly constrained at present. Observations aimed at measuring 
post-Newtonian effects would need to fit the local Galactic potential as well. Stars, however, 
are expected to be the dominant component of the extended galactic mass distribution 
near the central MBH. From observations, the stellar distribution around the MBH is best 
approximated by a double power-law density profile, with a (somewh at uncertain) break 



radius that may range from ~ 0.1 — 0.2 pc up to r ~ fewxlpc (ISchodel et al.l 12007 



20091 ). This translates into a single power-law model p(r) = po{r/ro)~'^ throughout the whole 
region of special interest for us, r < O.Olpc. Therefore, we adopt the following gravitational 
potential Vcai 



An 



(3-7)(2-7) 



.2-, M.(ro) /ry-^ 
Po^or ' = 7^ ( — ) , 7 7^ 2, 



(2 - 7)ro V'^o 



(38) 



In I — 



7 = 2, 



where M^^irn) is the total stellar mass within the radius tq. We will adopt tq = O.Olpc in 
this papero 

In the absence of definite observational measurements, and in order to chose the values 
for the slope 7 and for the normalization of the mass distribution, we have to appeal to 



^It is easy to see that n = 7 — 2, so that 1/2 < 7 < 3 imphes —3/2 < n < 1/2 in and (|37p. 
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stellar dynamical theory. The relaxation time in the Milky Way's nucleus is Tr ~ O(lGyr) 
and therefore old stellar populations may have had enough time to reach a relaxed steady- 
state. Mass segregation around a MBH leads to steep er profiles for heavy st ars {e.g. compact 
remnants such stellar black holes) than for light stars (IBahcall &: Wolilll977l). Recent F okker- 



Preto fc Amaro-Seoanell2009l ) show 



Planck and A^-body studies ( Alexander fc Hopman 20091 : 
that mass segregation is indeed a generic and robust property of the relaxed populations 
around a MBH; furthermore, it is stronger than expected according to Bahcall & Wolf, and 
leads to power-law density profiles with slopes such as 7heavy ~ 1.8 — 2.3 and 7iight ~ 1.0 — 1.6. 
According to the latter studies, the total amount of stellar mass packed inside ro = O.Olpc 
is M*(ro) ~ a X 2 X IO^Mq, where a = 0(1). Therefore, in our numerical tests, we will 
adopt 7 = 1.5,2.1 and a normalization for the total stellar mass M*(ro) = 2 x 1O^~^M0. 
This leads to circular velocities, induced by stellar mass alone, of order 20-50 km/s. Current 
observational constraints on the mass normalization are still roughly one or two orders of 



magnitude above of t hese values, whereas 7 is still very weakly constrained (IGhez et al.ll2008 
Gillessen et aDl2009h . 



Figures [3] and H] illustrate astrometric perturbations from the relativistic terms and 
from the local Galactic model (155]) . Not surprisingly, for an orbit with semi- major axis a 
the galaxy's perturbation effects are oc M^(a). For the adopted parameters, the cusp slope 
has only a weak effect but it is noticeable; furthermore, the galaxy's pertubations are also 
stronger for steeper cusps. Given the model adopted for the stellar cluster, the cumulative 
mass distribution is M(r) = M^,(ro)(r/ro)^~''', so when the slope 7 changes by an amount 
A7, the total mass within a (orbital semi-major axis) changes by 



AM (a) = M,(ro) 



> 3— 7— A7 

a \ fa 



To 



To 



3-7 



(39) 



Therefore, increasing 7 = 1.5 to 2.1 leads to a mass increment within a, AM*(< a) ~ 
0.19M*(a) ~ 3.8 x IO^Mq. Therefore a larger proper motion obtains which is, to first order, 
given by AZ2.1 ~ AZi,5 + A/i 5 ^ AM(a) ~ 0.89 mas. The agreement with the plots in the two 
bottom panels of Figure H] is very good. Figure H] also shows that the small difference (for 
an S2-like orbit) between the purely relativistic perturbation and the combined relativistic 
plus Galactic perturbation is smaller than astrometric capabilities — even in the case when 
the extended stellar cluster is relatively massive. 

Figure [S] shows the kinematic effect of relativistic and Galactic perturbations. As in the 
previous figures, the perturbations are shown as a function of mean anomaly. It can be seen 
that the relativistic kinematic perturbation on a S2-like star is strong enough to be detected 
by current instruments, for which 6v ~ 10 km/s. This is not the stronger case as the S14 
star has a closer pericenter passage and thus suffers a stronger (by a factor of ~ 5) kinematic 
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perturbation. 

In summary, three more things are evident from these figures that are worth commenting 

on. 

• For most of the known S-stars, the Galactic perturbations could be comparable to the 
relativistic effects, and even rather similar in time dependence. Perturbations on S2 
and Sf4 are, however, essentially dominated by relativistic effects. (Relativistic pertur- 
bations could dominate even more in stars further in, if such stars are discovered; this, 
combined with shorter periods, would make detecting the signature of the relativistic 
effects much neater.) It may be hoped, however, that even for the already-discovered 
stars, that the known and very specific form of the relativistic effects could enable 
them to be extracted, but it remains to be demonstrated. 

• Kinematics perturbations are concentrated near pericenter passage, whereas astromet- 
ric perturbations are of the same order throughout the orbit. For highly eccentric 
orbits, the the astrometric perturbation Al has two different origins: near pericen- 
ter, Al is due to Av, which shifts the phase of the orbit; near apocenter, Al comes 
precession Acu, amplified by the lever arm of the orbit. 

• Although the speed is maximal at pericenter, the kinematic perturbation is maximal 
(and can be much greater) at a phase before or after. 

5. Comparison with other post-Newtonian formulations 

In computing the preceding numerical results, we have made some non-trivial choices 
of convention. 

1. We have worked in what may be called be the Boyer-Lindquist gauge, which for zero 
spin reduces to the standard gauge of the Schwarzschild spacetime. 

2. Our independent variable is the affine parameter, which is proportional to the proper 
time. 

3. Our kinematic variables are the canonical momenta, and not derivatives of the coor- 
dinates. In particular, when computing Keplerian elements, we have fed momentum 
and not velocity values into the classical formulas. 

Different choices can lead to some surprising differences in the results. Of course, observable 
quantities must not change. But for abstract quantities such as the instantaneous a, even 
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the sign of the post-Newtonian effect can switch. We now explain how to convert to other 
conventions, considering for simplicity only the leading order Schwarzschild effects. 

For a weak-field Schwarzschild spacetime, the most common gauge choice is the harmonic 
form, where the metric is 

ds^ = ~ (l-- + X] dt^ + + (40) 
\ r J \ r J 

Working out \g^'^ we readily derive the Hamiltonian 

1/2 2 \ 1 / 2 



with e labeling orders as before, through the replacements 

r ^ e^^r p— i>ep. (42) 

Since there is no explicit dependence on t, the conjugate momentum pt will be constant. As 
before, the initial value of pt just sets the units of the affine parameter r, but if we set 

^ V 2 rj \2r^ 2r 8 J ^ ' 

initially then H = + O(e^), and the affine parameter will equal the proper time along the 
orbit, to the given order. 



As an aside, in Hamiltonian dynamics there is another possible interpretation of fl43|) : 
we can take the function pt(x, p) as a Hami ltonian in its own right, having three degrees of 
freedom and t as the independent variable. ISaha fc Tremaind (Il994j ) used this Hamiltonian 
to incorporate the leading-order post-Newtonian effects into a symplectic algorithm for long- 
term integration of planetary orbits. 



Considering now the coordinate velocity v = d'x./dt we have 

dr J ydr J \dpt J V^(^P 
which for the Hamiltonian (14T]) gives 



,,H#V7?^) = (^V (#t1 (44) 



v = j,rMi + -<'l 7i - -tM p + oCt"). (45) 



r 

Substituting from ( l43l) and rearranging, we have 

p2 3 



2 +-]e']p + Oie% (46) 
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Thus while x x p is a constan t of motion, in x x v the modulation P6l) appears (cf. the 
leading term in Equation 4.9 in lBlanchet fc Iyerll2003l . taking their u = in the test particle 
limit). 

As noted above, for the main results of this paper, we have computed orbital elements 
using the "momentum convention" , as 

-Ya-Y-l «(1-^)HxxpP (4T) 
and so on. On the other hand, if the "velocity convention" 

- — = — -- a(l-e^) = Ixx (48) 
2a 2 r ^ ^ ' ' ^ ^ 

is adopted instead, it follows from (146|) that the orbital elements can change at O(e^). 

Figure [6] shows the perturbations of the orbital elements (as in the top two panels of 
Figure [2]), but now computed by following the "velocity convention". It can be seen that the 
semimajor axis perturbation changes its sign; while u now shows a small oscillation around 
pericenter. 



6. Conclusions 

The prospect of detecting general relativistic effects in the S-stars near the Galactic 
center has recently aroused interest. To help gain more insight into these small but exciting 
effects, we have derived a simple formulation for the orbit equations and an algorithm for 
numerically integrating them, in which the post-Newtonian dynamics appears clearly as 
perturbations of the Kepler problem, and then examined the size of the perturbations. 

The post-Newtonian Hamiltonian is (fTOj) or equivalently (11511161) . These look like and 
are fairly simple generalizations of the Keplerian Hamiltonian, the main difference being 
that time is a coordinate and the affine parameter is an independent variable. For numerical 
integration we adapt the variable-timestep symplectic integrators recently developed for 
cometary or other highly eccentric orbits in the solar system. 

For photons a somewhat different approximation applies than for (comparatively) slow- 
moving stars, and we derive the Hamiltonian f|T2l) . It leads to gravitational redshift and 
related effects, but we leave the computation of these for future work. 

With the orbit integrator we proceed to compute the effects of the post-Newtonian terms 
on the orbital elements, sky position, and kinematics, taking the orbit of S2 as an example. 
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For S2, general relativity implies an astrometric effect of tenths of a mas and tens of km/sec 
in kinematics. Two surprising features are: (i) the astrometric effect is of the same order 
near pericenter as at apocenter, and (ii) the kinematic effect is greatest near (but not at) 
pericenter. 

Newtonian perturbations due to other masses in the Galactic-center region are unkown, 
but could plausibly be of the same order for the stars so far know. Disentangling the 
relativistic contribution from the Galactic perturbations may be the hardest problem in 
practice. It would be necessary to fit simultaneously for post-Newtonian effects and Galactic 
perturbations, without knowing the specific form of the Galactic perturbation in advance. 
Whether this is achievable is an open question. Simulations of the modeling pipeline are 
needed to get a clear answer, but the known and very specific time-dependence of the post- 
Newtonian effects suggest that we can be optimistic. 

We thank the referee, Chfford Will, for confronting some of our numerical results 
with linear perturbation theory, leading eventually to the comparison with standard post- 
Newtonian theory that forms Section 5 of the paper. 

MP acknowledges support by DLR (Deutsches Zentrum fiir Luft- und Raumfahrt). 
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Fig. 1. — Verification of the properties of the integrator. The upper panel shows the maxi- 
mum error in H over 2000 periods for four orbits with a = 2.4 x 10"^ and e = 0.5, 0.88, 0.95 
and 0.99. The lower panel shows the pericenter angle difference between the integration of 
an S2-like orbit, with initial time steps l/$(0) being 10~^ and 10~^ of the Keplerian period, 
and the theoretical value in (!29l) . The oscillations were averaged out from the time series 
to highlight its (linear) secular evolution, although a small oscillation remains in the former 
case. 
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Fig. 2. — Perturbation of the (osculating) orbital elements of an S2-like orbit (see Eq. [281) 
due to post-Newtonian terms. The upper panel shows a in gravitational units, the middle 
panel and lower panels show uo and Vt. The mean anomaly is taken as the independent 
variable, as a surrogate for time (see text). All angles are in degrees. 



- 22 - 



CO 



O 

CO 

E 
E 

Q) 

tn 



3 
o 

D) 
C 
10 

i 

CD 
-I— • 

c 

CD 

_o 

CD 
Q. 



CD 

•a 
o 
c 



24600 



24400 - 



24200 - 



24000 



-134.709 




-200 -150 -100 
-134.708 I r 



200 



-134.71 



-200 



-150 



-100 



-50 50 

mean anomaly M 



100 



150 



200 



Fig. 3. — Perturbation of the orbital elements of a S2-like orbit due to post-Newtonian 
terms plus a model for Galactic perturbations consisting of a stellar cluster around Sgr A* 
(see Eq. 1551) . The panels are analogous to the two upper panels in Figure [21 For the 
Galactic contribution, the larger the stellar-cluster and the steeper the cusp, the stronger 
the perturbation of the elements. The assumed total mass for the stellar cluster is, according 
to current theoretical estimates, relatively high. 
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Fig. 4. — Astrometric perturbations for an S2-like orbit. The upper panel shows Ab from 
post-Newtonian effects only, the middle panel from the model Galactic perturbations only, 
while the lower panel combines both perturbations. The precession due to both the ex- 
tended mass distribution and the relativistic effects is too small to be detected with current 
astrometric capabilities, even though the stellar cusp is relatively massive. 
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Fig. 5. — Kinematic perturbations. As in Figure HI the upper panel is from post-Newtonian 
effects only, the middle panel from the model Galactic perturbations only, and the lower 
panel combines both. Only the 10% of the orbit around pericenter passage is shown here. 
The kinematic perturbation due to PN terms on a S2-like orbit appears measurable with 
current spectroscopic resolution, bv ~ 10 km/s. 
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Fig. 6. — Perturbation of the (osculating) orbital elements of an S2-like orbit due to post- 
Newtonian terms as in Figure [21 but now computing the orbital elements using the "velocity 
convention" . There are noticeable differences with Figure O (i) the perturbation on the 
semi-major axis changes its sign; (ii) uj now shows a small oscillation around pericenter. 



